!************************************************************
! this subrouitne calculate invese wall point of isentropic
! supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-12
! Revised    :
!************************************************************
! Nomenclature: 
!   lm~V0   c-line properties, m-c-, p-c+, 0-streamline
!   lp~T02  coefficients of functions, *p-c+, *0(*)-streamline
!   x2f	   previous x2 value
!   x4f~V4f previous properties at point 4
!   temp1   interpolation coefficient of point 2, (x2-x1)/(x3-x1)
!   i,j     loop indicators
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine IWall()
   !         4()
   ! 3
   !    (2)
   !        1
   use VariableDef
   implicit none
   ! given variables
   ! real*8::x1,y1,V1,theta1,p1,rho1,x3,y3,V3,theta3,p3,rho3,x4,y4,theta4
   ! temporary variables
   real*8::lm,Vp,alphap,thetap,pp,rhop,ap,Mp,yp
   real*8::p00,rho0,V0
   real*8::lp,Qp,Sp,Tp,R0,A0,T01,T02
   real*8::x2f,p4f,rho4f,V4f
   real*8::temp1
   integer::i,j
   ! calculate variables
   ! real*8::V4,p4,rho4,T4,M4,Vv4,Vu4
   
   lm=(y3-y1)/(x3-x1)
   
   p4=(p1+p3)/2
   rho4=(rho1+rho3)/2
   V4=(V1+V3)/2
   p2=p4
   rho2=rho4
   V2=V4
   theta2=(theta1+theta3)/2
   
   p4f=0.0
   rho4f=0.0
   V4f=0.0
   
   do i=1,icor
      ! iteration for point 2
      x2f=0.0
      do j=1,icor
         Vp=(V2+V4)/2
         thetap=(theta2+theta4)/2
         pp=(p2+p4)/2
         rhop=(rho2+rho4)/2
         ! determine coefficients
         call Thermo(pp,rhop,Vp,ap,T4,Mp)
         alphap=asin(1/Mp)
         lp=tan(thetap+alphap)
         ! determine x2, y2 and theta2
         x2=(lm*x1-lp*x4-y1+y4)/(lm-lp)
         y2=y4+lp*(x2-x4)
         temp1=(x2-x1)/(x3-x1)
         theta2=theta1+(theta3-theta1)*temp1
         ! determine V2, p2, rho2 by interpolation
         V2=V1+(V3-V1)*temp1
         p2=p1+(p3-p1)*temp1
         rho2=rho1+(rho3-rho1)*temp1
         if (i==1) then
            p4=p2
            V4=V2
            rho4=rho2
         end if
         if (abs(x2-x2f) <= 1e-6) exit
         x2f=x2
      end do
      
      pp=(p2+p4)/2
      rhop=(rho2+rho4)/2
      Vp=(V2+V4)/2
      thetap=(theta2+theta4)/2
      yp=(y2+y4)/2
      ! determine coefficients along Mach line 24
      call Thermo(pp,rhop,Vp,ap,T4,Mp)
      alphap=asin(1/Mp)
      lp=tan(thetap+alphap)
      Qp=sqrt(Mp**2-1)/(rhop*Vp**2)
      Sp=delta*sin(thetap)/(yp*Mp*cos(thetap+alphap))
      Tp=-Sp*(x4-x2)+Qp*p2+theta2
      if (i==1) then
         p00=p3
         rho0=rho3
         V0=V3
      else
         p00=(p3+p4)/2
         rho0=(rho3+rho4)/2
         V0=(V3+V4)/2
      end if
      ! determine coefficients along streamline 34
      call Thermo(p00,rho0,V0,ap,T4,M4)
      R0=rho0*V0
      A0=ap**2
      T01=R0*V3+p3
      T02=p3-A0*rho3
      ! determine p4, V4, rho4
      p4=(Tp-theta4)/Qp
      V4=(-p4+T01)/R0
      rho4=(p4-T02)/A0
      ! justify if point 4 arrived
      if ((abs(p4-p4f) <= ep*p4f) .and. (abs(rho4-rho4f) <= erho*rho4f) .and.&
         &(abs(V4-V4f) <= eV*V4f)) exit
      
      p4f=p4
      rho4f=rho4
      V4f=V4
   end do
   call Thermo(p4,rho4,V4,ap,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)

end subroutine IWall
